^.o AFWL WUL.2) 

K/RTLAND AFB, N mex 


IMPULSIVE MOTION ON 
A FLAT PLATE PULSED 
WITH UNIFORM HEAT FLUX 

by Thomas H. Cochran and Eva T. Jun 

Lewis Research Center 
Cleveland, Ohio 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION • WASHINGTON, D. C. • MAY 1969 




TECH LIBRARY KAFB, NM 

IHNIH1 


0132035 


IMPULSIVE MOTION ON A FLAT PLATE PULSED 
WITH UNIFORM HEAT FLUX 

By Thomas H. Cochran and Eva T. Jun 

Lewis Research Center 
Cleveland, Ohio 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


For sale by the Clearinghouse for Federal Scientific and Technical Information 
Springfield, Virginia 22151 - CFSTI price $3.00 




ABSTRACT 


An analytical study was conducted for a semi -infinite flat plate that had fluid impul- 
sively moved on its surface and at the same time pulsed with uniform heat flux over a 
prescribed length. The numerical solution, which was obtained by an implicit finite- 
difference method with non uniform lattice dimensions, was compared with a solution for 
transient conduction in a slab and with a solution of the complete energy equation from 
the integral method of Karman-Pohlhausen. The solutions were confined to positions on 
the plate over which the flow is strictly transient and to fluids with Prandtl numbers 
greater than 1. 
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SUMMARY 

An analytical study was conducted for a semi-infinite flat plate that had fluid impul- 
sively moved on its surface and at the same time pulsed with uniform heat flux over a 
prescribed length. The numerical solution, which was obtained by an implicit finite- 
difference method with nonuniform lattice dimensions, was compared with a solution for 
transient conduction in a slab and with a solution of the complete energy equation from the 
integral method of Karman-Pohlhausen. The solutions were confined to positions on the 
plate over which the flow is strictly transient and to fluids with Prandtl numbers greater 
than 1. 


INTRODUCTION 

In recent years, transient forced-convection heat -transfer has become of increasing 
interest to scientists and engineers. Solutions for this class of problems have applica- 
tions to devices, such as the rocket engine and the nuclear reactor, in which startup and 
shutdown are important phases in the operation cycle of the equipment. The need for 
solutions to transient forced-convection problems may also be necessitated by experi- 
mental testing in facilities that afford relatively short testing times, such as the shock 
tunnel and the zero-gravity drop tower. In the former, "pulse type" testing is conducted 
while in the latter, experiments may not be initiated until after they are released. 

Transient forced-convection heat-transfer problems may be classified as external 
(flat plate) or internal (pipe) flows. A further breakdown denotes the transient as hydro- 
dynamic, thermal, or both hydrodynamic and thermal. In the hydrodynamic problem, 
thermal boundary conditions are not a function of time and the fluid flow is unsteady. 
Transient thermal problems are characterized by constant fluid flow, and the wall heat 
flux or the wall temperature are a prescribed function of time. The third case, of 
course, treats both unsteady flow and unsteady thermal boundary conditions. A survey 
of the work in this field has been provided by Soliman (ref. 1). Recent work includes 
references 2 and 3. 

The purpose of this report is to solve an external problem that is unsteady both hy- 
drodynamically and thermally. The problem may be described as follows: A flat plate 
is initially in thermal equilibrium with quiescent fluid that is above it. At some pre- 


scribed time, the fluid is impulsively set in motion, and the plate receives a step in- 
crease in heat flux over a portion of its surface. Rozenshtok in a recent paper (ref. 4) 
treated a similar problem except that he assumed a step increase in surface tempera- 
ture rather than in heat flux. He also did not consider the case for an unheated entrance 
length. In the present work, the problem is solved analytically by use of the integral 
method of Karman-Pohlhausen. These solutions are compared with a numerical solution 
which is obtained by use of an implicit finite-difference method with nonuniform lattice 
dimensions and which is assumed to best approximate true conditions. The validity of 
this assumption is substantiated by the close agreement between the numerical solution 
and an analytical solution for transient conduction in a slab for times and positions on the 
plate for which conduction is the heat-transfer mechanism. The solution is confined to 
positions on the plate over which the flow is strictly transient and to fluids with Prandtl 
numbers greater than 1. 


ANALYSIS 

The equations to be solved for transient laminar flow in a thermal boundary layer 
are 

Continuity: 


^ + 9Y = 0 

ax ay 


Momentum: 


SU+U— +V — - _ ii 

at ax 3y ay 2 [i dx 

Energy: 

0T + n 9T t y 9T _ v /a 2 T\ f v /3U\ 2 

at ax ay Pr^ 0y 2 J c p \a y/ 


( 1 ) 


( 2 ) 


(3) 


where the symbols are defined in appendix A. The fluid is assumed to be incompressible 
and the viscosity constant. Because of these restrictions, the continuity and momentum 
equations are independent of the energy equation and, therefore, may be solved directly 
for the velocity components U and V. 
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Fluid Flow Conditions 


The fluid flow problem associated with impulsive motion on a flat plate has been 
treated by numerous investigators, references 5 and 6 being representative. Their solu- 
tions require simplification of the governing equations and are of two types. One result 
is based on a linearization of the momentum equation to an Oseen-type (Rayleigh’s meth- 
od) and another utilizes the momentum -integral method (Karman-Pohlhausen solution). 
Both solutions indicate that there are two domains on the plate that are not analytically 
joined. In the first domain, the velocity profile is solely a function of time; in the sec- 
ond domain, the velocity profile is only dependent on distance from the leading edge of 
the plate. A primary difference between the two results is the prediction of the extent of 
the domains. The Rayleigh solution indicates that it occurs at x = U^t (ref. 5), whereas 
the momentum -integral result shows that it occurs at x- 0.385 U^t (ref. 6). All symbols 
are defined in appendix A. 

With a knowledge of the aforementioned work, Stewartson (ref. 7) treated the full 
time -dependent momentum equations and obtained solutions in the form of power series. 
His results indicate that for x > U^t the time -dependent Rayleigh solution is the correct 
solution and that for x < 0. 377 U^t the position-dependent momentum -integral solution 
is the appropriate form to use. The domain between these two, U^t > x > 0. 377 U^t, is 
a transition region in which the velocity is dependent on both time and position. These 
results are shown graphically in figure 1. A detailed discussion of Stewartson's work is 
contained in reference 8. 



Domain 

I 

Transient 

II 

Transition 

III 

Steady state 

x* Uoot 



Figure 1. - Hydrodynamic conditions for impulsive motion on flat plate. 


3 



The present work is confined to those positions on the plate for which x > U^t. 
Therefore, an exact solution for the velocity profile is available and is 


V = 0 


U = U^erf (- 1 — \ 

\2 yjvtj 


(4) 

(5) 


However, for the purpose of obtaining an analytical solution for the energy equation, the 
velocity should be in a form that can be conveniently integrated. An expression in ref- 
erence 4 that closely approximates equation (5) is the fourth-degree polynomial 


U = 





( 6 ) 


where = 3. A comparison between equations (5) and (6) is made in figure 2. 



Figure 2. - Approximate and exact transient velocity profiles for impulsive motion on 
semi-infinite flat plate. 
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Thermal Conditions 


The plan to be followed in finding a solution for the thermal problem is to approach 
it from both an analytical and a numerical viewpoint. The primary analytical approach 
used is the integral method whereas an implicit finite -difference method was used to find 
the numerical solution. A comparison of the solutions from the different methods can 
then be made. 

Integral solution . - The thermal problem to be solved is that of a semi-infinite flat 
plate subjected to liquid impulsively set in motion in the plane of the plate that at the same 
time is pulsed with a uniform heat flux over a portion of its surface (i.e. , the plate has 
an unheated entrance length x Q ). Attention is confined to those positions in the heated 
region for which the initial hydrodynamic transition x = U TO t is upstream. 

The energy equation (eq. (3)), when viscous dissipation is neglected and equation (4) 
is applied, becomes 


0T n 9T v 9 2 T 

9t 3x Pr gy2 

The boundary conditions are 

T = T^ at (x Q , y, t) 

T = T m at (x > x Q , °°,t) 

— = - — at (x > x , 0, t) 

3y Ak ° 

and the initial condition is 


(7) 


T = T„ 


at (x > x Q , y, 0) 


The problem is solvable by the integral method. Therefore, equation (7) becomes 
■5m = /*6 

1 U(T - T Jdy = - 

Pr 9 y l y =0 


— Z' 6 ' 1 ’ (T - T )dy + A f &T U(T - T 0 

at y o ax J o 


,)dy = — — — 


( 8 ) 


where 6 T is defined as that position above the surface of the plate where T = T^ and 
Pr is the Prandtl number /iC p /k. Equation (8) is written for fluids with Prandtl num- 
bers greater than 1, and the solution will be confined to this case. 
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The formulation may be further simplified by defining a new dependent variable 


e = 


T - T 


OO 



The boundary and initial conditions become 

6 = 0 at (x Q , y, t), (x>x Q , 5 T ,t), and (x > x Q , y, 0) 

— = -1 at (x >x 0,t) 

3y 

and equation (8) is simplified to 

— f T 6 dy + JL /* 6t XJ9 dy = ~ (9 

at •'0 3x J 0 Pr 

The dependent variable is now approximated by a polynomial in y. An expression 
that has been used by others working on transient thermal problems (see refs. 9 and 10) 
and has yielded acceptable results is 



( 10 ) 


Substitution of equations (6) and (10) into equation (9) yields 


! 96 t 2 

6 at 



v 

Pr 


(H) 


where E = 6rp/6jj. For the particular case at hand, E < 1 so that the term in paren- 
thesis may be approximated by an arithmetic average with little error: 


1 9 ( 6 t 2 ) + 13 



3(6 t 3) 


6 at 


168 


ax 


v_ 

Pr 


(12) 


This equation may be simplified by letting 


6 t 2 = 


so that 


1 El aw + -i! u °° Pr m 1 / 2 w = i 

6 v 3t 112 Si 


(13) 


Equation (13) is solvable by the method of characteristics, the details of which are pre- 
sented in reference 11. The system of equations for the characteristics of (13) may be 
written as 


6 



The solutions of the system (14) are 


6 


T 



6 t = 2.04 


v 6 H^ “ x o } 


U TO Pr 


1/3 


x-x Q = 


0.468 


(Pr) 


1/2 


Uoot 


(14) 


(15) 


(16) 


( 17 ) 


Substitution of equations (15) and (16) into equation (10) results, respectively, in 



r 


e = 1.02 


y 6 H<* - x o } 


^ > 11/3 


U^r 


1 - 


v. 


0. 490 y 

s 

y 6 h (x - x Q ) 

1/3 

U TO Pr 



(19) 


The solution predicts two regimes on the plate that are not analytically joined. The 
temperature in one (eq. (18)) is dependent only on time (conduction regime) while in the 
other (eq. (19)) it is dependent on distance from the thermal leading edge, time, and free- 
stream velocity (convection regime). 

The analysis indicates that the boundary between the conduction and convection re- 
gimes is described by equation (17). However, in light of Stewartson's hydrodynamic 
analyses, it seems reasonable to assume that there is a transition regime between these 
two regimes. For this particular solution, equation (17) describes the upper extent of the 


Regime 

A Conduction 
B Transition 

u oo C Convection 


x = 0.377 u^t x=uj. x = F 2 u 00 t + x x^FiU^t + x,. 



Figure 3. - Thermal conditions for impulsive motion on flat plate with 
unheated entrance length. 


convection regime whereas a conservative estimate of the lower limit of the conduction 
regime is x - x Q = U^t. A better estimate for this limit is determined later in the sec- 
tion Extent of Regimes. The conditions existing on the plate are shown graphically in fig- 
ure 3 where the extent of the regimes is denoted in a general form. 

Since the solution obtained by this method depends on the form of the polynomial that 
was assumed for 9 in terms of y (see eq. (10)), other solutions to the problem are pos- 
sible by assuming different polynomials. A form in common use besides the parabolic 
profile is the quartic profile (ref. 4). For the problem at hand, the boundary conditions 
necessary to form the required equation are 
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o 

II 

CD 

at y 

- 6 t 

— = o 
ay 

89 = 

at 

y = 0 

d 2 e ( 

ay 



ay 2 " 

2 

at 

y = o T 


3y 





= 0 at y = 0 


The quartic polynomial is 


9 = 



( 20 ) 


and the solutions to the problem are 


6 = 1.290 



- y + 0. 150 


(?) 



( 21 ) 


e = 1 . 10 


- V 

1/3 

- y + 0o 207 

U^Pr 

2/3 

v 2 - 0 047 

U.Pr 

U TO Pr 


6 h ^ - x Q )_ 

J \J m % 

v 5 H 0f - x Q )_ 


( 22 ) 


x-x 0 = 


0.443 


Uoot 


[(Pr) 1 / 2 

Numerical solution. - The form of the energy equation to be solved numerically is 


(23) 


— + U M erf 71 = a — (24) 

at ax g y 2 

where rj = y/2 a is the thermal diffusivity and x - x - x Q . An approximate solution 
of equation (24), with its initial and boundary conditions, was found using an implicit 
finite -difference method. The partial differential equation was replaced by a finite- 
difference equation which was then solved. This process is equivalent to finding 0(x, y, t) 
at discrete lattice points in the x-, y-, t-space (fig. 4). Because of the exponential nature 
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t 



Figure 4. - Lattice dimensions in x, y, t-space. 


of the solution, and our interest near the origin, nonuniform lattice dimensions were used. 
These dimensions may be expressed as 


Ax. = x. - Xj_ x - PXINTCAx.^) 

Ay k = y k -y k -i = PYINT ( A yk-i) 


where PTINT, PXINT, and PYINT are greater than or equal to one. 

The numerical values of PTINT, PXINT, PYINT, At 2 , Ax 2 , Ay 2 , KMAX, and J MAX 
were determined by trial and error under the condition that the solution not be affected by 
lattice sizes, that 8 rea t er than or equal to the heated length of the plate, and 

that y^MAX £ rea t er than or equal to the maximum estimated thermal layer thickness. 

The finite -difference equivalent of equation (24), as derived in appendix B, is 


0 j,k-l,n ( ~ p ) 


+ e 




A' 


+ y + p( 1 + e)J + 6 


j, k+1, n 


(-Pe) = W 


j, k, n-1 + y0 j-l,k,n 


(25) 


where e, j3, y, and p are defined in appendix A, and for simplicity 0(x.,y, , t ) is ex- 

] K n 

pressed as 6 . , . Initial and boundary conditions are 

h n 


0 j,k, 1 = 0 over a11 i» k 
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over all k, n 


0 l,k,n-° 


0 j,KMAX+l,n 


= 0 overall j,n 


30 

3y 


y=o 


Ay. 


( 0 , 


j,2,n j, l,n 


di 1 J = - 1 j > 1, n > 1 


If equation (25) is written for all lattice points of y > 0 at fixed values of x and t 


t 



Figure 5. - One row of lattice points. 


(see fig. 5), a system of KMAX - 1 number of linear equations with KMAX - 1 unknowns 
is obtained: 



1 


where Z. = 0. . - , m = KM AX - 1, and a, b, c, and d are defined in appendix B. 

This system of equations is a tridiagonal matrix, all elements of which are zero except 
for the main diagonal and the two diagonals on either side of it. The coefficients of Z, 
the values for which may be found in appendix B, are functions of lattice point positions. 

At the first assumed time t = t 2 and displacement x = Xg, the elements on the right side 
of the system (26) are determined by the initial and boundary conditions. The variable 0 
is now determined from the matrix by Gauss' elimination (see ref. 12). The solution of 
the matrix is presented in appendix C. With 9 known for the first x-row of lattice points, 
another set of KMAX - 1 equations can be written for the next row of lattice points in the 
same time plane. The elements on the right side of the matrix for this solution are ob- 
tained from the results of the solution from the preceeding x-row. This process is con- 
tinued until the maximum displacement in x is reached. After obtaining the solution for 

all lattice points in the x, y-plane at the first time, 0. . can now be found at successive 

J 5 n 

values of time by repeating the above processes until the time reaches the limit of inter- 
est, TMAX. 

A FORTRAN IV computer program was written to find the solution of the finite- 
difference equation in the manner just described. A description of the computer program 
with flow charts and a listing of the programs is given in appendix D. 

Transient conduction solution . - The integral solution indicates that for some posi- 
tions on the plate the mode of heat transfer is strictly conduction. This suggests that for 
such positions a solution may be obtained from the energy equation in the form 


3T 

at 


" v \ a 2 T 
Pr ' ay 2 


with the boundary conditions 


T = T^ at (y, 0) 


(27) 


T = T^ at (°°, t) 


3T _ %/_ 

dy Ak 


at (0, t) 


This problem has been solved in reference 13 and the solution is 
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RESULTS 


Conduction and Convection Regimes 

A comparison between the analytical and numerical solutions was made for water as 
a working liquid. The conditions assumed were a free-stream temperature of 98. 9° C 
and free-stream velocities of 3.05 and 12.2 centimeters per second. The surface heat 
flux was taken as 1261 watts per square meter. In comparing the results, the numerical 
solution was assumed to always best approximate the true conditions. 

The validity of the conduction equations was tested by assuming velocities, times, and 
positions such that x - x Q > U^t. The results are shown in figures 6(a) and (b). The 
curves indicate that, of the integral solutions, the parabolic profile is in best overall 



(a) Free-stream velocity in x-direction, 3. 05 centimeters per second; distance from leading thermal edge, 
0.915 centimeter; time, 0.154 second. 


Figure 6. - Temperature profiles. 


3.5 



50x10 
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(b) Free-stream velocity in x-direction, 12.2 centimeters per second; distance from 
leading thermal edge, 2. 72 centimeters; time, 0.116 second. 



(c) Free-stream velocity in x-direction, 3.05 centimeters per second; distance from leading thermal edge, 
0.915 centimeter; time, 2.02 seconds. 


Figure 6. - Continued. 




(d) Free-stream velocity in x-directlon, 12.2 centimeters per second; distance from leading thermal 
edge, 2. 72 centimeters; time, 1.52 seconds. 


Figure 6. - Concluded. 

agreement with the numerical solution, the greatest discrepancy between the two occur- 
ring at large distances from the surface. This is understandable in light of the finite 
boundary condition imposed on the integral solution at large y as opposed to the infinite 
one that actually exists. The close agreement between the transient conduction solution 
and the numerical solution is expected and validates the numerical approach to the prob- 
lem. 

The convection equations were evaluated for conditions so that x - x Q < 0. 324 U^t 
(eq. (23)). This criterion was chosen from the two available ones (eqs. (17) and (23)) be- 
cause it minimized the size of the convection regime, thus ensuring that the selected con- 
ditions were definitely out of the transition regime. The results in figures 6(c) and (d) 
show that the quartic profile is in the best overall agreement with the numerical solution. 

In summary, for the conduction regime, the transient conduction solution (eq. (28)) 
is representative of the temperature profile. For the convection regime, the integral so- 
lution using a quartic profile (eq. (22)) provides a good approximation for the temperature 
distribution. 


Transition Regime 

The nature of the temperature distribution between the conduction and convection re- 
gimes was investigated by selecting a velocity (3. 05 cm/sec) and a position (0. 915 cm) 
and determining the numerical solution for different times. The results are presented in 
figure 7 where the conduction and convection solutions, as determined in the previous 
section, are also plotted. 
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The conditions in figure 7(a) are such that the position selected is just out of the con 
duction regime, as defined by the extent criteria 00 > x - x Q > U^t. As seen from the 
figure, the conduction solution is still generally in good agreement with the numerical 
solution. However, at large distances from the surface, the two begin to deviate from 
one another, signaling the start of convection effects. This result may be logically ob- 
tained by considering the physical processes that are occurring. Conduction dominates 
as the mode of heat transfer at a position in the thermal boundary layer until fluid that 
was originally upstream of the thermal leading edge reaches the position in question. 
Since the liquid velocity is greater at large distances from the surface, the convection ef 



(a) Free-stream velocity in x-direction, 3.05 centimeters per second; distance from leading thermal 
edge, 0. 915 centimeter, time, 0.643 second. 



(b) Free-stream velocity in x-direction, 3.05 centimeters per second; distance from leading thermal 
edge, 0.915 centimeter; time, 0.856 second. 


Figure 7. - Temperature profiles. 
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fects are realized first at those distances. 

The results in figure 7(b) are for a later time than those in figure 7(a). Here again, 
the conduction solution is in fairly good agreement with the numerical solution; however, 
the deviation between these two solutions at large distances from the plate is greater than 
that in figure 7(a). Comparison of the convection curve in figures 7(a) and (b) indicates, 
that as time progresses, it comes into closer agreement with the numerical solution. It 
is obvious then, that the importance of convection as a mode of heat transfer increases 
with time. 
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(c) Free-stream velocity in x-direction, 3.05 centimeters per second; distance from leading thermal edge, 
0.915 centimeter; time, 1.516 second. 


Figure 7. - Concluded. 


For figure 7(c), the conditions are such that the selected position is just in the con- 
vection regime as determined by the extent criterion 0<x-x Q <0.324 U^t. The cur- 
ves indicate that the convection solution is now in best agreement with the numerical solu- 
tion and that the transition from conduction to convection as the mode of heat transfer 
has penetrated to small distances from the surface of the plate. 

Two primary conclusions are to be drawn from these results: The first is that during 
this regime the mode of heat transfer is changed from conduction to convection. The 
second conclusion is that in the transition regime the conduction solution of equation (28) 
is a good approximation for the temperature distribution. 
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Extent of Regimes 


From the analysis, it was determined that the upper extent of the convection regime 
was defined by equation (23). The function Fg in the extent criterion, as shown in fig- 
ure 3, may then be expressed as Fg = 0. 443/(Pr)*/ 2 . For the conduction regime, its 
lower extent has been conservatively estimated as x - x Q = U^t. As discussed in the 
previous section, the onset of convection effects at a position (x - x Q ) in the thermal 
boundary layer is signaled by the arrival at the position of liquid that was originally up- 
stream of the leading thermal edge. Since the velocity boundary layer is larger than the 
thermal boundary layer for fluids with Prandtl numbers greater than 1, the liquid in 
transit from the leading thermal edge to a position is not continuously moving at the free- 
stream velocity. Therefore, x - x Q = U^t is conservative in that it overestimates the 
size of the transition regime. 

A criterion for the lower extent of the conduction regime may be obtained by consid- 
ering liquid particles in laminar flow at some distance y above the surface of the plate. 
An expression for the rectilinear motion of these particles is 


x 



U dt 


(29) 


where x m denotes the displacement over which the liquid particles moved at the free- 
stream velocity and t^ the time required to move this distance. Inserting equation (6) 
into equation (29) and integrating results in 


X 0 = X oo + U oo 


' 1. 096 y\ / t 1/2 

V ,1/2 A 



0.082 y 3 \ / 1 
,3/2 )[ t l/2 



f 0. 0056 y 4 V 1 

v , 2 N 



The displacement x^ may be expressed as 

X oo = Uootoo 


(30) 


(31) 


where t^ is determined from the time required for the hydrodynamic boundary layer to 
increase in size to the distance y under consideration. A relation for the hydrodynamic 
boundary layer, as presented in the section Fluid Flow Conditions, is 
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(32) 


6 h = 3.65(^t) 1 / 2 

so that 

0. 0752 Uy 2 
= (33) 

v 


The initial convection effect at a position (x - x Q ) is realized when liquid at a distance y 
reaches the position at the same time that the conduction thermal boundary layer has in- 
creased in size to y. If a criterion for the thermal boundary layer is assumed to be 


T - T„ 


= 0.01 


T W " T ® 


equation (28) may be solved such that 
Inserting equations (32) to (34) into equation (30) results in 


; \l/2 


(34) 


x ' x o 


3. 507 2.687 _ 0.587 _ 4. 608 

.(Pr) 1 / 2 (Pr) 3 / 2 (Pr) 2 Pr 


(35) 


where F ^ (see fig. 3) is the term in parentheses. A plot of the extent criteria of equa- 
tions (23) and (35) is provided in figure 8. 



Figure 8. - Extent of regimes as function of Prandtl number. 
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SUMMARY OF RESULTS 


An analytical and numerical study of the thermal conditions on a flat plate with an 
unheated entrance length, subjected to fluid impulsively moved, and pulsed with uniform 
heat flux yielded the following results: 

1. For positions on the plate denoted as the conduction regime and describable by the 
conditions x > U^t and F jU^t + x < x < °°, where x is the distance measured from 
the leading edge along the plate, is the free stream velocity, t is time, Fj is a 
function of the Prandtl number that is characteristic of the lower limit of the conduction 
regime, and x Q is the unheated entrance length, the temperature profile may be pre- 
sented as 


T - T 


= 2Q W /j*\ 
Ak \Pr/ 


1/2 


ierfc 


2 f—\ 

[PrJ 


1/2 


where 


3.507 f 2.687 _ 0.587 _ 4.608 
(Pr) 1 / 2 (p r ) 3 / 2 Pr 2 Pr 


for Pr > 1 


and T is the fluid temperature, T M is the free stream temperature, is the surface 
heat flux, A is the surface area, k is the fluid thermal conductivity, v is the kinematic 
viscosity, Pr is the Prandtl number, and y is the distance measured perpendicular to 
the plate. 

2. For positions on the plate denoted as the transition regime and describable by the 
conditions x > U^t and FgU^t + x Q <x < FjU^t + x Q , where Fg is a function of the 
Prandtl number that is characteristic of the upper limit of the convection regime, the 
temperature profile may be approximated as 


2 Qw f vt\ ^/ 2 

T - T = — - (—1 ierfc 
Ak \Pry 



1/2 


where 
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for Pr > 1 


F _ 0.443 

2 " (Pr) 1 ^ 

3. For positions on the plate denoted as the convection regime and describable by the 
conditions x > U^t and x Q < x < + x Q , the temperature profile may be best pre- 

sented as 


V (x ~ x o } 
U_Pr 


y + 0.207 


UooPr 

" x o) 


y - 0.047 


y 5 h (x - x 0 ) 


where 5 jj is the hydrodynamic boundary layer thickness. 

Lewis Research Center, 

National Aeronautics and Space Administration, 
Cleveland, Ohio, February 27, 1969, 
124-09-17-01-22. 
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APPENDIX A 


SYMBOLS 


A 

a,b, c, d 

S 

E 
F i 


k 

m 

P 

Pr 

PXINT 

PYINT 

PTINT 

Q 

q 

T 

TMAX 

t 

U 

V 

XMAX 

x 

x 

YMAX 

y 


2 

heated area of plate, cm 

array of constants in a tridiagonal matrix 

specific heat, J/(kg)(°C) 

ratio of thermal to hydrodynamic boundary layer, 

function of Prandtl number that is characteristic of lower limit of conduction 
regime 

function of Prandtl number that is characteristic of upper limit of convection 
regime 

thermal conductivity, J/(m)(sec)(°C) 

order of matrix equation, KMAX - 1 
o 

pressure, N/m 
Prandtl number, pC /k 

r 

growth rate of x intervals, Ax. /Ax. ^ 

growth rate of y intervals, Ay^Ay^j 

growth rate of t intervals, At n /At n _j 

heat flux, W 

total velocity, cm/sec 

temperature, °C 

maximum time of interest 

time, sec 

velocity in x-direction, cm/sec 
velocity in y-direction, cm/sec 
total heated length, cm 

distance measured from leading thermal edge along plate, x = x - x Q , cm 
distance measured from leading edge along plate, cm 
maximum estimated thermal boundary layer thickness, cm 
distance measured perpendicular to plate, cm 
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Z unknown in matrix equation representing 0 at all vertical points above plate at 

fixed t and x 

o 

a thermal diffusivity, cm /sec 

P 1/At n 

y erf rj/Ax. 

6 boimdary layer thickness, cm 


e 

*7 

0 

M 

y 

P 

* 


1/PYINT 

argument of the error function, 
dependent variable in energy equation, 0 = 
dynamic viscosity, cP 

9 

kinematic viscosity, (m) /sec 

a/(Ay k ) 2 

5 t 2 


T - 

Qvy/Ak 


Subscripts: 

H hydrodynamic conditions 

JMAX number of lattice points on x- scale 

j,k, n indices denoting specific lattice points on the x-, y-, and t-scales, respectively 
KMAX number of lattice points on y-scale 

o thermal entrance condition 

T thermal conditions 

W conditions at y = 0 

°° free-stream conditions 
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APPENDIX B 


FINITE-DIFFERENCE APPROXIMATION 


The finite -difference approximation of the energy equation in the form 


— + U^erf (r?) ^ = a 
dt dx 



is derived in the following manner: From figure 4, it may be seen that 


de 

at 


j,k,n 


At ^j,k,n “ 0 j,k,n-l) 


80 

dx 


\h K n 


k n 

Ax. ] ’ K ’ n 


0 j-l,k, n> 


The second partial derivative in y is obtained by considering two consecutive intervals 


h-\ 



on the y- scale (see fig. 9): 


a 2 g 

av 2 




(0. , - 29. , +9. , , ) 

v J,k+€,n j,k,n J,k-l,n' 


where 9. , can be expressed by a linear interpolation, 9. . = eQ. , , 

3 , K+€ ? n 9,9 j,K+t,n 

+ (1 - e)0. i_ . Although the existence of 3 6/ 0y proves that 0. t cannot be ac- 
j ? k, n K+t ? n 

curately approximated by a linear interpolation, for the case when e = 0. 91, the results 
derived from the linear interpolation differs from that of a second-order approximation 
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by only 0.2 percent. The partial differential equation in finite -difference form is then 


x U^erffa) 

~ t ~ ( 0 j,k,n " 0 j,k,n-P + ~ ^j,k,n " 0 j-l,k,n^ 


? _ e0 j,k, +1, n ~ ^ + e ) 0 j,k,n + 0 j,k-l,n] 




with initial and boundary conditions 


V,i = ° 


0 l,k,n- 0 


0 j, KMAX+1, n “ 0 


Ay<- 


(0 j,2,n ” 0 j,l,n 


j > l,k> 1 
k > l,n > 1 
j > i, n > i 

)=-l j > 1, n > 1 


Equation (Bl) can be rearranged in the form 


e 


i,k-l,n<-P) + 9 J,k,n[^ ^ + p(1+ £) ] + e J,k + l,n ( - p€ > = pe j,k,n-l + r Vl,k,n 


(B2) 


where p,/3,y, and e are defined in appendix A. Applying the boundary conditions at 
k = 2 yields 

0 j,k,n^ +y + pe) + 0 j,k+l,n (_pe) = /30 j,k,n-l + y0 j-l,k,n ~ p(0 j,k,n ~ 0 j,k-l,n) (B2a) 

The last term on the right side is p^ygM-l). Applying the boundary condition at 
k = KM AX gives 


e j,k-l,n<-' > ) + f ’j,k,„[/ ! + r + P(1 + e >] ^ f 'i,k,n-l + 


(B2b) 


For every given Xj and t n , equations (B2) make up a system of KMAX - 1 equations 

with KMAX - 1 unknowns, 6. , (k = 2, 3, . . .KMAX), in the following form: 

J> k j n 
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b^Zj + c^Zg “ d l 

a 2 Z l + b 2 Z 2 + C 2 Z 3 _ d 2 

a 3 Z 2 + b 3 Z 3 + C 3 Z 4 = d 3 


a m Z m-l + b m Z m ' d m 


where 


Z i ^j, i+1, n 


m = KM AX - 1 


bj = /3 + r + pe 


C 1 - c 2 c m _i - -P e 

d l = ^j,2,n-l + ^j-l,2,n " P^YaX' 1 ) 


a 2 - a 3 - . . . - a m - -p 


b 2 = b 3 = • • • = b m = ^ + y + p(l + e) 


d i (i=2 ’ 3 ’ ’ • * ' m > = &j, i+1, n-1 + ^j-1, i+l, n 
0. 1 is obt a ined by a pplying the boundary condition 

J? *9 11 

^j,2,n " 9 ), 1,t) = " 1 


Ay.- 


where 


<, i,l,n = 5 i ,2,n- a y2<- 1 > 


(B3) 
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APPENDIX C 


SOLUTION OF TRIDIAGONAL MATRIX EQUATION 

The system of m tridiagonal equations (eq. (26)) can be solved by using Gauss' 
Elimination to yield the results 

Z = d* 
m m 


Z. = d? - c*Z. , 
1 1 ii+I 


i = m - 1, m - 2, . . .2, 1 


where 


Ci 

r * - 1 
C i - — 
b l 


d 1 

d l=- 

b l 


A 


* 

C i = 


C i 


b i " a i C i-l 


> i- 


d * d i- a i d *-l 

i_b i-vr-i 


i = 2, 3, . . . m 




and c T and d*(i = 1, 2 , ... m) are computed first and then used to find Z^(i = m, m - 1, 

... 1 ). 

If (bj - a j c *_j) = 0, c? and d* are not computed but rather 


d. - a.d? * 

j* l l l-l 

d i+l = 7 


c r+i = ° 


z i+i = d r+i 


and 


Z i = 


d i+l ~ b i+l Z i+l ~ c i+l Z i+2 
a i+l 
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APPENDIX D 


COMPUTER PROGRAM 

The following sections give a detailed description of the computer program written to 
solve equation (24) numerically. The program consists of a main program and a sub- 
routine FINDZ. Included herein are the general description of the program, a dictionary 
of the FORTRAN symbols used, a flow diagram (fig. 10), and the program listings. The 
input/output format is described, and the computer input/output of a sample case is also 
presented. 


Program Description 

Main program . - The input data are read, and the computation is begun by initializing 
x, y, t and inserting the boundary and initial conditions. For each subsequent row in x, 
the distance along x is incremented, the constants and coefficients of its matrix equa- 
tion are calculated, and subroutine FINDZ is called. 

Subroutine FINDZ . - The root of the matrix equation Z is found by the method de- 
scribed in appendix C . The column matrix Z is returned to the main program through 
blank COMMON. 


Dictionary of FORTRAN Symbols 


FORTRAN Engineering 

symbol symbol 


X 

Y 

T 

XMAX 
JMAX 
KM AX 
DELX 
DELY 
DELT 


X " X o 


y 

t 


A ^k 

At n 


distance along plate -x 

array containing lattice dimensions in y -direction 
time 

maximum x of interest 

number of lattice points in x-direction 

number of lattice points in y-direction 

j_i_ 

length of J in increment in x 
length of K increment in y 

in 

length of N in increment in t 
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FORTRAN 

Engineering 


symbol 

symbol 


DELXO 


value of 1. /PXINT times the first x-increment 

DELYO 


value of l./PYINT times the first y-increment 

DELTO 


value of l./PTINT times the first t-increment 

THETA(J, K, 2) 

k, n 

evaluated at interval in time 

THETA(J, K, 1) 

0 j,k,n-l 

at the N in interval in time 

PARTYO 


partial derivative of 0 with respect to y at y = 0 

RNU 

V 

kinematic viscosity, input constant 

ALPHA 

a 

thermal diffusivity, input constant 

BETA 

P 

i./Atn 

GAMMA 

y 

U^erf n/ Ax. 

SQ 


2 yjvt 

CD2 


(a/Ayg) times PARTYO 

C2 


variable equal to 0 or CD2, used in computing the 
D -array 

UINF 

u ro 

free-stream velocity, input constant 

EPS 

e 

i./pyint 

ERFARG 

V 

y/2-^t 

FF 


array containing values of U^erf^) at all y-intervals 

M 

m 

KM AX - 1 

A, B, C,D 

a,b, c, d 

arrays containing coefficients and constants of eq. (26) 

cs(D 

«? 

defined in appendix C 

DS(I) 

d i 

defined in appendix C 

DENOM 


B(I) - (Aft) times CS(I - 1)) 

ITAG(I) 


I in element of array ITAG, an indicator of whether 
DENOM is zero; 0-no, and 1-yes 

Z 


array output from subroutine FINDZ, 
Z(K) = THETA(J, K+l, 2) 
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KLAST 


on exit from FINDZ, KLAST is last nonzero element of 
Z ; KLAST is immediately set equal to KLAST + 1 
denoting the last nonzero element of THETA(J, K, 2) 
at given J 

JTV, NTV because of large volume of results, solution at every 

lattice point was not written out; solutions were 
written out only at every (NTV) th interval in time 
and (JTV)^ interval in x 

result of MOD(N-l, NTV) 

result of MOD(j-l, JTV) 

result of 1. -THETA(J-1,KLAST,2)/THETA(J,KLAST,2); 
at small values of time, when x becomes greater 
than some specific X(J), dependent on value of T, 
THETA no longer changes significantly with x; when 
RAT is sufficiently small, THETA(JN, K, 2) is set 
equal to THETA(J, K, 2) for all values of K and for 
all J < JN < JMAX 

name list name containing values of ALPHA, RNU, 
U1NF, PARTYO, XMAX, TMAX, JMAX, KM AX, 
DELXO, DELYO, DELTO, PX1NT, PYENT, NTV, 
JTV 

Blank COMMON contains M, A, B, C, D, Z, KLAST 
The program listing is as follows: 

C **** IMPLICIT FINITE DIFFERENCE EQUATION *********************** 

C *************** IRREGULAR LATTICE SIZES ************************************ 

C RESULTS WRITTEN OUT ONLY AT EVERY NTV INTERVALS OF T AND EVERY JTV 

C INTERVALS OF X 

COMMON M, A,B,C,D,Z,KLAST 

DIMENSION A(IOO), B(IOO), C(IOO), DI100), ZtlOO), YUOO)* THETA I 12 
11, 100,2), FF(IOO) 

NAMELIST/ INPUT/ AL PHA , RNU, U I NF , PART YO , XMAX , T MAX , JM AX , KMAX , DELXO , 

1DELY0, DELTO, PXINT,PYINT,PTI NT, NTV, JTV 
LOO READI5, INPUT) 

M=KM AX - 1 

Y ( 1 ) = 0. 

CD2= -ALPHA * PARTYO/(DELYO*PYINT) 

EPS = l./PYINT 
N= 1 
T = 0. 

DELT = DELTO 
DO 3 K= 1 , KMAX 
DO 3 J= 1 , JMAX 
3 THET A ( J , K , 2 ) = 0. 

C AT T=0. , THETA=0. FOR ALL X AND Y 

20 N •= N + 1 

DELT = DELT * PTINT 
T = T + DELT 

IF (T .GT. TMAX ) GO TO 100 


MN 

MJ 

RAT 


INPUT 


31 



n o 


00 10 J= l , JM AX 
00 10 K= 1 jKMAX 

10 THET A( J, K * 1 ) = THETAIJ»K,2) 

C THETA(J,K,2) STANDS FOR THE TA I J , K, N ) , THETA(J,K,1> FOR THETA ( J,K ,N— 1 1 

MN = MOD C N— 1 »NTV ) 

IF ( N .NE. INTV+1) .AND. MN .£Q. 0) WRITEI6,64) T 

BETA = l./OELT 

SQ = 2.* SORT (RNU*T ) 

X * 0. 

OE LX = DELXO 
DO 12 K= l »KMAX 

12 THETA! 1 , K , 2 ) = 0. 

C AT X * 0., THE TA=0. FOR ALL Y AND T 

DO 11 J = 2 * JMAX 
MJ = MODI J-l » JTV ) 

DELX = PXINT * DELX 
X = X + DELX 

IF (X .GT. XMAX) GO TO 20 

X AND T ARE NOW FIXED, COMPUTE COEFFICIENTS AND CONSTANTS 
OF FINITE-OIFFERENCE EQUATION AT EACH POINT OF Y 
DELY = DELYO 
DO 1 K = l , M 

IF (J .GT. 2) GO TO 2 
IF (N .GT. 2) GO TO A 
C DELY,. A, C, AND Y ARE INDEPENDENT OF X AND T, COMPUTED ONLY ONCE IN 

THE PROGRAM, AT J=2, N=2 

C FF I K ) IS INDEPENDENT OF X, COMPUTED ONLY WHEN J = 2 AT EVERY VALUE OF N 

DELY = PYINT * DELY 
ROE= ALPHA/OEL Y**2 
YIK+l) = Y ( K ) + DELY 
AIK) = -ROE 
CIK) = -ROE*EPS 
t* ERFARG = YIK + D/SQ 

FF I K ) = ERF I ERF ARG ) *UINF 
2 GAMMA = FF I K ) /DELX 

BIK) = BETA + GAMMA - I A I K ) +C IK) ) 

D2 = 0. 

IF (K-l) 6,6,1 

6 D2 = CD2 

BIK) = BETA + GAMMA - CIK) 

AIK) = 0. 

1 DIK) = BETA * THETA I J * K + l > 1 ) +GAMMA*THETAt J-l , K + l, 2) + 02 
CIM) = 0. 

IF (J .EQ. 2 .AND. N .EQ. 2) WRITE(6,66) I Y < K ) ,K= 1 , KMAX ) 

IF I N .EQ. INTV+1) .AND. J .EQ. (JTV+l) ) WRITEI6.63) T 
CALL FINDZ 

C KLAST IS THE LAST NONZERO ELEMENT OF THE SOLUTION FROM FINDZ 

KLAST = KLAST+1 
DO 5 K=1 , M 

5 THETA! J,K+1,2) = ZIK) 

THET A I J , 1 , 2 ) = THETA(J»2»2) - I P Y I N T*DE L YO ) * PART YO 

IF I MN .EQ. 0 .AND. MJ .EQ. 0 ) HR I TE I 6 , 65 ) X , ( THE T A I J , K , 2 l , K= 1 , KL AST ) 

13 RAT = 1. - THETA! J-l, KLAST, 2)/ THETA! J, KLAST, 2) 

IF (ABSIRAT) .LT. l.E-04) GO TO 7 

C IF THETA NO LONGER CHANGE WITH X, JUMP OUT OF DO 11 COOP 

11 CONTINUE 
GO TO 20 

7 JP = J + 1 

DO 8 JN = JP , JMA X 
DO 8 K=1 , KMAX 

8 THETAIJN, K, 2) = THETA ( J ,K,2> 

IF I MN .EQ. 0 ) WRITE(6,67) X 
GO TO 20 
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>£) ® -g 0 s U1 


' ■-y"3SK®;^ 

63 FORMAT ( IHL,55X, 9HAT TIME = # F8.5, 2X, 4HSEC • ) 

64 FORMAT ( !Hlt55Xf 9HAT TIME = , F8.5, 2X, 4HSEC.) 

65 FORMAT ( / ,30X, 6HAT X = , E13.5, 3H CM, 2X, 59HTHETA AT VARIO 

LUS VERTICAL POINTS ( NONUN I FORML Y D I STR LBUTED ) , / < 1C < IP6 l 3. 4 ) ) > 

66 FORMAT < LH 1 * 45X, 4 1H VALUES OF Y AT WHICH THETA ARE CALCULATED# / 

1 ( 10{ 1PE13.4) ) ) 

67 FORMAT ( 1HL # 30X, 48HTHET A NO LONGER CHANGES WITH X AT X GREATER TH 
1 AN , E 1 3 • 5 , 3H CM) 

END 


$I8FTC TRID 

SUBROUTINE FINDZ 

C SOLUTION OF TRIDIAGONAL MATRIX EQUATION BY GAUSS^ ELIMINATION 

C REFERENCE — FORSYTHE AND WASOW P.104 

D I MENS I UN A( 10 0) ,B( LOO ) , C ( 1 00 ) , D ( 1 00 ) , Z ( 1 00 ) , C S ( 10 C ) • D S ( 100 ) , 
1ITAGC 100) 

COMMON M, A,B,C,D,Z,KLAST 

ITAGf 1 ) = 0 

CS ( 1 ) = C ( 1 >/B ( 1 ) 

DS(1) = D(l) / B ( 1) 

DO 4 I =2 » M 
I T AG ( I ) -=0 

IF ( I T AG ( 1-1 ) .EQ. 1 ) GO TO 4 
DENOM = BID - A ( I )*CS< 1-1) 

IF (ABS(DENOM) - l.E-37) 1,1,3 

1 ITAGf I ) = 1 

C DENOM I I ) = 0. NO CS(I) AND NO DStI) COMPUTED 

DS(I*1) = (DII)-A(I)*DSU-1))/C(I) 

CS ( 1+1) = o. 

GO TO 4 

3 DS ( I ) = l D ( I ) - A(I)*DS(I-1) ) /DENOM 
CS!!) = C ( I ) /DENOM 

4 CONTINUE 
DO 5 1=1, M 
K=M-I+1 

IF ( AB S ( D S ( K ) ) .GT. l.E-37) GO TO 6 
Z I K ) = 0. 

Z(K) = DS ( K ) 

KLAST = K 

KLAST IS THE LAST NON-ZERO Z ELEMENT 
KM 1 => K-L 
DO 9 J = 1 , KM 1 
I = K— J 

IF (ITAG(I)) 7,7,8 
Z( I ) = DS ( I ) - CS< I )*Z< 1+1) 

GO TO 9 

Z ( I ) = f D ( I + 1 ) — B( I + 1 ) * Z ( I + l) — C ( I + 1 ) *Z ( I +2 ) ) / A ( I + 1 ) 

CONTINUE 
RETURN 
END 
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Description of Input and Output 


Input. - Input data are read into the program in namelist form. A namelist statement 
in the main program specifies a list of variables under the namelist name INPUT. Data 
cards must be written in the form appropriate for namelist input for the computer used. 
For the present case, the computer used is the IBM 7090/7040 direct-couple system, and 
the input card format is specified in reference 14. The first card of each case must have 
a $ sign in the second column followed immediately by the namelist name. Each input 
variable is given in the following form: variable name = number, separated by commas. 

A $ sign must follow the last data item indicating the end of a namelist data group. In the 
program, the system of units used for input/output is not restricted. In the input data 
cards for the following sample case, all the values of the dimensioned physical variables 
are given in cgs-centigrade units: ALPHA and RNU in square centimeters per second, 
UINF in centimeters per second, XMAX, DELXO, and DELYO in centimeters, and TMAX 
and DELTO in seconds. All other input variables are dimensionless quantities. The fol- 
lowing listing shows an input data card of a sample case: 


SINPUT AL PHA = 1 .67226—03 , RNU=3 . 1 586 E-03 , UINF=3.048, PARTY0=-1.0, 

XMAX=5.0, TMAX= 1 . E— 04 « JPAX=121, KMAX=100, DELXC=5. 5E-09, DEL Y0*4 . 5E-06 , 
DELT0= 1. E-05» PX I NT= 1 . 2 * PYINT=1.1, PTINT = 1.1, NTV=3, JTV=12$ 


Output . - The lattice positions in the y-direction are written out at the beginning of 
the computation. Values of THETA at these y-positions will be written out at increasing 
values of x at every NTV^ 1 value of t. 

In some instances, one may wish to see the results in terms of T - T^ instead of 
6 defined presently as (T - T oc )/(Q w /Ak). Rather than multiplying the final result at each 
point by the constant (Q^/Ak), it is more expedient to redefine 6 as T - T^. The differ- 
ential equation will remain the same. The only alteration required in the program is the 
input 39/ 3y at y = O(PARTYO), which, with the new definition of 6, has to be multiplied 
by ((^/Ak). The following listing shows the computer output of a sample case: 
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values of y at which theta are calculated 

o. 4.9 500E-06 1.0395E-05 1.6304E-O5 2.2973E-05 3.0220E-05 3.8192E-05 4.696IE-05 5.6608E-05 6.7218E-05 
7 • 8890E-05 9.1729E-05 1.0585E-04 1.2139E-04 1.38486-04 1.57276-04 1.7795E-04 2.0070E-04 2.2572E-04 2.5324E-04 
2-8351E-04 3.1681E-04 3.5344E-04 3.9374E-04 4.3806E-04 4.8682E-04 5.4045E-04 5.9944E-04 6.6434E-04 7.3572E-04 
8. 1425E-04 9.0062E-04 9.9563E-04 1.1001E-03 1.2151E-03 1.3416E-03 1.4807E-03 1.6337E-03 1.8020E-03 1.9872E-03 
2.1908E-03 2.4149E-03 2.6613E-03 2.9324E-03 3.2306E-03 3.5586E-03 3.9194E-03 4.3163E-03 4.7529E-03 5.2331E-03 
5*76 13E-03 6« 3424E-03 6.9816E-03 7.6847E-03 8.4582E-03 9.3089E-03 1. 02456-02 1.1274E-02 1.2407E-02 1. 36526-02 
1.50Z2E-02 1.6530E-02 1.8107E-O2 2.00L1E-02 2.2 017E-02 2.4224E-02 2.6651E-02 2.932LE-02 3.2258E-0 2 3.5489E-02 
3.9043E-02 4.2952E-02 4.7252E-02 5.1983E-02 5.7106E-O2 6.2909E-02 6.9205E-02 7.6131E-02 8.3749E-02 9.2128E-02 
1.01 35E-01 1.1 1496-0 l 1.2264E-01 1.3491E-01 1.4840E-01 1.6325E-01 1.79586-01 1.9754E-01 2.1730E-0L 2.3904E-01 
2.6294E-01 2.8924E-01 3.1B17E-01 3.5000E-01 3.8500E-01 4.2351E-01 4.6586E-01 5..12456-01 5.6370E-01 6.2008E-01 


AT TIME = 0.00004 SEC. 


AT X * 0.26L23E— 06 CM THETA AT VARIOUS VERTICAL POINTS ( NONUN IFORMLY DISTRIBUTED) 


6.7824E-05 

6.2874E-05 

5.7472 E- 05 

5. 16336-05 

4 ■ 5398E-05 

3.0852E-O5 

3.21 34E—05 

2.5446E-05 

1.90616-05 

l ■ 3298E-05 

8.4738E-06 

4.81 64E-06 

2.3762E-06 

9.0811E-O7 

3.3636E-07 

9.1274E-O0 

1 .9312E-08 

3.1335E-09 

3.8546E-10 

3 .57016-1 1 

2. 4816E-12 

1.2941E-13 

5.071 6E- 15 

1.4985E-16 

3.3533E-18 

5.7155E-20 

7 .46776—22 

7 • 534 IE— 24 

5.9156E-26 

3 • 6453E-2 8 

U7783E-30 

6.9267E-33 

2.1 715E-35 










AT X = 0, 

, 25904E-05 CM 

THETA AT VARIOUS VERTICAL 

POINTS (NONUNIFORMLY DISTRIBUTED) 


1.42 18E-04 

1.3723E-04 

1. 3181 E- 04 

1.2588E-04 

1. 1943E-04 

1.1242E-04 

l .0484E-04 

9.6685E-05 

8. 7985E-05 

7.8785E-05 

6. 91 75E-05 

5.9292E-05 

4.9336E-05 

3.9569E-05 

3.0316E-05 

2.1940E-05 

1.4789E-05 

9.1294E-06 

5 . 05926—06 

2.461 IE-06 

1.0263E-06 

3.5837E-07 

1 • 0256E-07 

2. 3633E-08 

4.3260E-09 

6.2 3636-10 

7.O409E-11 

6.2425E— 12 

4.3412E-13 

2.3816E-14 

1-0369E-15 

3.6083E-17 

1 .0108E-18 

2.2953E-20 

4. 2498E-22 

6.4466E-24 

8.0379E-26 

8.25406-28 

6 . 9909E — 30 

4.8051 E— 32 

2.8168E-34 

1.3216E-36 











AT X = 0. 

. 2 3357E-04 CM 

THETA AT VARIOUS VERTICAL 

POINTS (NONUNIFORMLY DISTRI 

1 BUTED ) 


2.4571 E-04 

2.4076E-04 

2. 3536 E— 04 

2.2950E-04 

2.2313E-04 

2.162 3E-04 

2 - 08 76E-04 

2.00706-04 

1.9201E-04 

1.8269E-04 

l. 7271E-04 

1 . 6209E— 04 

1.5082E-04 

1. 3895E-04 

1.2653E-04 

1.1364E-04 

1 . 004 IE— 04 

8. 69966-05 

7. 36246-05 

6.055BE-05 

4.81 L8E-05 

3.6655E-05 

2.6524E-05 

1.8031E-05 

1. 1367E-05 

6.5478E-06 

3 • 3926E-06 

1.5 555E-06 

6.2141 E— 0 7 

2. 1335E-07 

6. 2254E-08 

1.531 3E-08 

3.1 5 B5E-09 

5.4454E— 10 

7. 83446- 1 1 

9.3967E-12 

9. 3881 E- 13 

7.8051E-14 

5. 3939E-15 

3.0949E-16 

1.4720E-17 

5.8O04E-19 

l . B969E— 20 

5.1271E-22 

l. 1464E-23 

2. 1200E-25 

3.24126-27 

4. 09616-29 

4.2784E-31 

3. 6929E-33 

2.6321E-35 

1.5045E-37 










AT X = 0.20852E-03 CM THETA AT VARIOUS VERTICAL POINTS ( NONUN LFORMLY DISTRIBUTED) 


2.7582E-04 

2.7087E-04 

2.6550E-04 

2.59686-04 

2. 5338E-04 

2 .46586—04 

2.3926E-04 

2.3 1 38E-04 

2 • 2295E-04 

2- 1394E-04 

2- 04 36E-04 

1.9420E-04 

1.8349E-04 

1.7225E-04 

1.60546-04 

1 .48406 — 04 

1 .35946-04 

1.2 326E-04 

1. I047E-04 

9 • 7 73 86-05 

8.5220E-05 

7. 3097E-05 

6.155 2E— 05 

5.07646-05 

4.08976-05 

3. 2089E-05 

2.4437E-05 

1.7992E-05 

1 . 2748E-05 

8.6459E-06 

5.5791 E— 06 

3.40106-06 

L. 94296-06 

L .03096-06 

5 • 0303E-07 

2.2341E-07 

8.9338E-08 

3. 1816E-08 

9.98 35E-09 

2.7319E-09 

6. 4575E-10 

1.3072E-10 

2.2492E-1 l 

3.2678E-12 

3. 9873E- 1 3 

4. 06736-14 

3.45586-15 

2.4384E-16 

1 • 4255E- 1 7 

6 • 8925E- 1 9 

2. 7525E-20 

9.069 IE — 22 

2.4636E-23 

5.5 1 50E-25 

1.017OE-26 

1.5446E-28 

1 .9319E-30 

1 .98996-32 

1.6875E-34 

l. 1031E-36 



AT X * 0. 

. 1 B594E — 02 CM 

THETA AT VARIOUS VERTICAL 

POINTS (NONUNIFORMLY DISTRIBUTED) 


2. 7593E-04 

2.7 G98E-04 

2.6560E-04 

2 . 59 7 BE -04 

2.5349E-04 

2.46696-04 

2.3936E-04 

2.3149E-04 

2. 2306E-04 

2. 1406E-04 

2. 0447E-04 

1.9432E-04 

1 . 836 IE- 04 

1. 72 38E-04 

1. 6067E-04 

1 .4054E-O4 

1 .36 096-04 

1.23426-04 

1 . 1065E-04 

9. 792BE-05 

8 ■ 5429E— 05 

7.3327E-05 

6. 1 806E-05 

5. 1044E-05 

4. 1206E-05 

3.24266-05 

2 .4 800E- 05 

I .83 74E-05 

1.3I40E-05 

9.O336E-06 

5. 9450E-06 

3.7275E-06 

2.215 3E-06 

1 .24126 — 06 

6.51696-07 

3.1867E-07 

1 .4417E-07 

5 . 9 92 1 E-08 

2.27I7E-08 

7.797 BE-09 

2. 4052E-09 

6.61586-10 

1.6107E-10 

3.44586-1 1 

6.4318E-12 

1.04056-12 

l .4498E-13 

1.73Q0E-14 

1 . 7587E-15 

1 . 5161E- 16 

1. 1037E-17 

6.7603E-19 

3.4732E-20 

1.4926E-21 

5. 3526E-23 

l. 59866-24 

3 .96936-26 

8. 18166-28 

I .398 26-29 

1.9792E-31 

2. 3183E-33 

2.2 4516-35 

1.7677E-37 










THETA NO 

LONGER CHANGES 

WITH X AT X 

GREATER THAN 

0. 555236-02 

CM 




AT TIME = 0 .00008 SEC. 

AT X = 0.26123E-06 CM THETA AT VARIOUS VERTICAL POINTS l NONUN I FORMLY DISTRIBUTED) 

7.0279E-O5 7.3329E-05 6.7915E-05 6.20356-05 5.570BE-05 4.8984E-05 4.1954E-05 3.4767E-05 2.7633E-05 2.0829E-05 

L. 4677E-05 9.4959E-06 5.5188E-06 2.8087E-06 I.2174E-06 4.3672E-07 I.2623E-07 2.8721E-O0 5.0486E-09 6.7627E-10 

6.8393E-U 5.1929E-L2 2.9521E-13 1.2558E-14 4.00L9E-16 9.5745E-18 1.7254E-19 2. 3514E-21 2. 4349E-23 1.9262E-25 

I.L712E-27 5.5096E-30 2.0201E-32 5.8153E-35 I.3040E-37 

AT X = 0.25904E-05 CM THETA AT VARIOUS VERTICAL POINTS I NONUN I FORMLY DISTRIBUTED) 

1.6655E-04 1.616 0£— 04 1.5617E-04 1.5022E-04 I.4371E-04 1.3660E-04 1.2888E-04 1.2051E-04 I.IL5LE-04 1-0I67E-04 

9. 1658E-05 8.0952E-05 6.9895E-0S 5.8692E-05 4.7617E-05 3.70I8E-05 2.7289E-05 1.8829E-05 1.1968E-05 6.B760E-06 

3.4935E-06 1.5326E-06 5.6652E-07 1.7235E-07 4.2267E-0B 8.2162E-09 1.2503E-09 1.4772E-10 1.3493E-11 9.5165E-13 

5.19O0E-14 2. 1 9766— 15 7.2591E-17 1.08306-lB 3.8642E-20 6.3218E-22 0.3O95E-24 8.8390E-26 7.6573E-28 5.4302E-30 

3. 1638E-32 1.5176E-34 5.9151E-37 


AT X = 0.23357E-04 CM THETA AT VARIOUS VERTICAL POINTS ( NONUN 1 FORMLY DISTRIBUTED) 


3.2403E-O4 

3.1 908E-04 

3. 1365E-04 

3 .07 7 IE-04 

3.0121E-04 

2.9409E-04 

2.8632E-04 

2.7703E-O4 

2.6859E-04 

2 .50546-04 

2.4763E-04 

2 . 358 IE-04 

2. 23076-04 

2 .09376-04 

1.94 726-04 

1.7915E-04 

1 .6271E-04 

1.45546-04 

1.2 70OE— 04 

1.0975E-04 

9. 17 406-05 

7.4202E-05 

5. 7648E-05 

4.2632E-05 

2. 9682E-05 

1 .9201E-05 

1.13866-05 

6.0504E-06 

2.8438E-06 

1. 1583E-06 

4.O175E-07 

I. 1686E-07 

2.01716-08 

5 . 58 1 8E-09 

9.04846-10 

1 .1982E-10 

1 .29716-1 1 

1 • 1 504E- 1 2 

8.3818E-14 

5 . 02996- 1 5 

2.4909E— 16 

1.0192E-17 

3.447 IE— 19 

9.6308E-21 

2.22806-22 

4.2566E-24 

6. 7207 E— 26 

8 . 78 85E-2 8 

9.45286-30 

8.4 198E-32 

6. 1962E-34 

3.7 L86E-36 











AT X * 0. 

20852E-03 CM 

THETA AT VARIOUS VERTICAL 

POINTS INONUN IFORMLY DISTRIBUTED) 


4.2695E-04 

4 . 2 200E-04 

4. 1660E-04 

4.10 7 IE-04 

4.0430E-04 

3.9732E-04 

3 . 8973E-04 

3.BI50E-04 

3. 72596-04 

3.6294E-04 

3. 5253E-04 

3.41 32E-04 

3 .292 7E— 04 

3- 1637E-04 

3.02 59 E — 04 

2 . 8794E— 04 

2. 7241E-04 

2 .56046—04 

2. 3889E-04 

2.2102E-04 

2-02 54E-04 

1.8 359E-04 

1 • 6435E— 04 

1.4504E-04 

l. 2590E-04 

1.0722E-04 

B.9316E-05 

7.2502E-05 

5.71006-05 

4.3399E-05 

3. 1632E-05 

2. 1942E-05 

1.4360E-05 

8.7 770E-06 

4 . 95 56 E -06 

2.5545E— 06 

1.1879E-06 

4.92626-07 

U8016E-07 

5. 7522E-O0 

1 .58 85 6-08 

3.76246-09 

7.5852E-10 

1.2930E-10 

1.8529E-U 

2. 2214E-12 

2.21916-13 

1.841 3E -14 

1 . 2656E- 1 5 

7. 19236-17 

3.3740E-18 
3.9390E— 36 

1 ■ 305 IE — 19 

4. 1588E-21 

1.0911E-22 

2. 3560E-24 

4. 1857E-26 

6. 1 180E-28 

7. 3567E-30 

7.2777E-32 

5.92336-34 



AT X = 0. 

18594E-02 CM 

THETA AT VARIOUS VERTICAL 

POINTS (NONUNIFORMLY DISTRIBUTED) 


4. 2860E-04 

4.2365E-04 

4 • 1825E-04 

4. 1236E-04 

4.0595E-04 

3.9897E-04 

3 -9140E-04 

3.83 18E-04 

3.74276-04 

3.6464E-04 

3.5426E-04 

3.4307E-04 

3.3106E-04 

3.1 820E-04 

3.0447E-04 

2.8 98 8E-04 

2 .744 36—04 

2 .58 16E-04 

2.4112E-04 

2.2339E-04 

2 ■ 0507E-04 

1.063 IE-04 

1 .67296—04 

l • 4822E-04 

1.29356-04 

1.10956-04 

9. 3331 E— 05 

7.6705E-O5 

6. 1600E-05 

4.8025E-05 

3. 6249E-05 

2.6376E-05 

1.04 1 5E— 05 

1 .22 726-05 

7 . 7622E-06 

4 .63106 — 06 

2.5887E-06 

1 . 3461 E — 06 

6.4627E-07 

2. 8420E-07 

1. 1354E-07 

4. 08646-08 

1.31 39E-08 

3-7426E-09 

9. 3682 E- 10 

2-0446E-10 

3.062OE-11 

6.2704E-12 

8 . 696 IE- 1 3 

1.0243E-13 

1.0196E-14 

6.5373E-16 

5.9904E-17 

3.5103E-18 

1.7129E-19 

6 .94326 — 2 1 

2.3331E-22 

6.4876E-24 

1.49O0E-25 

2. 8279E-27 

4.4239E-29 

5. 7032E-3 1 

6 . 0559E— 33 

5.28916-35 

3.6797E-37 








THETA NO 

LONGER CHANGES 

WITH X AT X 

GREATER THAN 

0.66628E-02 

CM 
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